function Sol = compute_eq_values(Par,vphi,mz,agg,mlabor,vztildeH,vzH_expost)
% Computes equilibrium values and store them in a structure
% Output:
%   Sol is the structure that stores all equilibrium values
% Input:
%   Par - parameter structure
%   vphi - 1 by ngoods vector of parameter phi for each Home firm
%   mz - 2 by ngoods productivity matrix
%   agg - 1 by 3 vector of solution for the aggregate variables
%   mlabor - 4 by ngoods matrix of labor allocation
%   vztildeH - 1 by ngoods vector of cutoff productivity
%   vzH_expost - 1 by ngoods vector of realized productivity shocks for Home
    ngoods = Par.ngoods;
    nRBgoods = Par.nRBgoods;
    nSBgoods = Par.nSBgoods;
    sigma = Par.sigma;
    rho = Par.rho;
    tau = Par.tau;
    kappa = Par.kappa;
    sigmazH = Par.sigmazH;
    zbarH = log(mz(1,:))-sigmazH^2/2;
    % Record wage rate in Foreign
    Sol.wF = agg(3);
    % Vector of dummies, = 1 when there is strike
    Sol.vstrike = (vzH_expost<vztildeH);
    % Actual rate of strikes
    Sol.strikerate_RB = sum(Sol.vstrike(1:nRBgoods))/nRBgoods;
    Sol.strikerate_SB = sum(Sol.vstrike(nRBgoods+1:end))/nSBgoods;
    Sol.strikerate_overall = sum(Sol.vstrike)/ngoods;
    % Theoretical probability of strikes
    Sol.strikeprob_RB = mean(logncdf(vztildeH(1:nRBgoods),zbarH(1:nRBgoods),sigmazH));
    Sol.strikeprob_SB = mean(logncdf(vztildeH(nRBgoods+1:end),zbarH(nRBgoods+1:end),sigmazH));
    Sol.strikeprob_overall = mean(logncdf(vztildeH,zbarH,sigmazH));
    % 4 by ngoods matrix of labor allocation at firm-market level
    Sol.mlabor = mlabor;
    % Ex post output at firm-market level
    Sol.moutput = mlabor.*kron([vzH_expost.*(1-kappa).^Sol.vstrike;mz(2,:)],[1;1])./[1;tau;tau;1];
    Sol.sector_output = [(Sol.moutput(1,:).^((rho-1)/rho)+Sol.moutput(3,:).^((rho-1)/rho)).^(rho/(rho-1));(Sol.moutput(2,:).^((rho-1)/rho)+Sol.moutput(4,:).^((rho-1)/rho)).^(rho/(rho-1))];
    Sol.final_output = sum(Sol.sector_output.^((sigma-1)/sigma),2).^(sigma/(sigma-1));
    Sol.price_index = ([agg(1);agg(2)]./Sol.final_output).^(1/sigma);
    Sol.sector_price = ([agg(1);agg(2)]./Sol.sector_output).^(1/sigma);
    % 4 by ngoods matrix of firm-market level freight on board prices
    Sol.mprice = repmat(Sol.sector_price.*Sol.sector_output.^(1/rho),2,1)./Sol.moutput.^(1/rho)./[1;tau;tau;1];
    % 4 by ngoods matrix of market shares
    Sol.mmarketshare = ([1;tau;tau;1].*Sol.mprice./repmat(Sol.sector_price,2,1)).^(1-rho);
    % 1 by ngoods vector of Home firms' profits before paying rents to the unions
    Sol.profitsH = sum([1;tau].*Sol.mprice(1:2,:).*Sol.moutput(1:2,:) - (1-kappa).^(Sol.vstrike).*mlabor(1:2,:),1);
    % 1 by ngoods vector of Foreign firms' profits
    Sol.profitsF = sum([tau;1].*Sol.mprice(3:4,:).*Sol.moutput(3:4,:) - agg(3).*mlabor(3:4,:),1);
    % 1 by ngoods vector of requests made by Home unions
    Sol.request = compute_request(Par,vphi,[vztildeH;mz(2,:)],agg,mlabor);
    % 1 by ngoods vector of rents paid to the unions
    Sol.rent = (1-Sol.vstrike).*Sol.request + Sol.vstrike.*vphi.*Sol.profitsH;
    % 2 by ngoods matrix of investment
    Sol.minvestment = [Sol.profitsH - Sol.rent;Sol.profitsF].*Par.s;
end